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Atmospheric  Structure  Simulation:  An  ARMA 
Model  for  Smooth  Isotropic  Two-Dimensional 

Geophysical  Power  Spectra 


1.  INTRODUCTION 

Atmospheric  fluctuations  in  wind  speed,  temperature,  and  density  are  characterized  by 
continuous  power  spectral  density  functions.  For  example,  two  dimensional  wind  speed  PSD's  are 
found  to  have  log-log  slopes  of  about  -2.67.  Such  spectra  often  are  used  In  simulating  an  envi¬ 
ronment  or  predicting  atmospheric  structure.  Fast  Fourier  transform  analysis  provides  a  means 
for  filtering  white  noise  with  spatial  filters  to  simulate  a  stationary  time  or  spatial  data  set  In 
many  applications  the  Fourier  transform  technique  provides  adequate  processing  speed.  At¬ 
mospheric  structure,  which  Is  multidimensional  and  extends  over  large  volumes,  is  not  readily 
simulated  by  conventional  techniques.  For  example,  stationary  baseline  Fourier  simulation  of  a 
two-dimensional  vertical  sheet  using  numerical  fast  complex  Fourier  transform  (FFT)  algorithms 
requires  -  4MNLog2(MN)  multiplications  on  the  two-dimensional  random  number  set  The 
resulting  Fourier  transform  values  then  must  be  multiplied  by  the  desired  power  spectrum  given 
the  required  vertical  (L^)  and  horizontal  (L^)  correlation  values  and  variances  (o2)  corresponding 
to  a  particular  altitude.  One  must  then  compute  the  inverse  2-D  transform  of  this  array.  The  total 
number  of  computer  multiplications  for  each  altitude  is  0{4MN(Log2(MN)+lll.  Since  L^,,  L^,  and 
o  vary  with  altitude,  M  2-D  sheets  must  be  calculated;  resulting  In  a  combined  number  of  mul¬ 
tiplications  of  0{4M2NlLog2{MN)+ 1 J) .  If  one  wishes  to  simulate  a  detailed  square  vertical  plane 
environment  having  a  side  of  100  km  with  100  m  resolution,  then  M>N>  1024  and  the  number 
of  multiplications  exceed  9.0  x  1010.  A  three  dimensional  analog  of  this  procedure  would  require 
(X4M4(3Log2(M)+l))  computer  multiplications  or  0(1.4  x  1014).  Even  high  speed  computers 
available  today  (-100  Mflop/sec)  would  take  -  378  hours  or  2.25  weeks  to  produce  a  single 
realization.  Thus  alternative  multidimensional  techniques  available  through  modem  spectral 
estimation  (and  perhaps  new  developments  In  chaos  or  wavelet  theory)  must  be  exploited  to 
reduce  the  computational  burden. 
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The  Phillips  Laboratory  Strategic  High  Altitude  Atmospheric  Radiance  Code  (SHARC) 1  uses 
first  principles  to  calculate  point  to  space  and  limb  viewing  atmospheric  background  Infrared 
radiance  and  transmittance.  Real  atmospheric  Infrared  background  perturbations  occur  from 
fluctuations  In  temperature  and  density  of  the  contributing  molecular  species.  Version  4  of  the 
SHARC  code  envisions  a  capability  to  evaluate  radiance  structure  from  estimated  variances  in  the 
standard  temperature  and  density  profiles.  This  report  studies  the  possibility  of  producing  two- 
dimensional  synthetic  structure  from  autoregressive/moving  average  (ARMA)  analysis  as 
contrasted  with  the  Fourier  method.  The  analysis  Is  performed  with  a  view  toward  producing 
circularly  symmetric  power  spectral  densities  that  account  for  Isotropic  correlations  in  the 
horizontal  plane,  including  the  line-of-slght.  that  existing  methods  approximate. 

During  the  course  of  the  study,  it  became  apparent  that  low-order  auto- regressive  (AR) 
analysis  alone  was  unable  to  produce  two-dimensional  power  spectra  for  the  large  correlation 
lengths,  evident  in  the  upper  atmosphere,  that  matched  desired  geophysical  specifications.  The 
study,  therefore,  was  extended  to  Include  full  ARMA  analysis.  Although  low  order  AR  analysis 
could  efficiently  be  used  near  ground  level  where  correlation  lengths  typically  are  small,  the  larger 
high  altitude  correlation  lengths  called  for  complete  ARMA  analysis,  where  structure  simulations 
required  application  of  a  high  order  recursion  technique.  As  part  of  an  ongoing  study,  this  report 
describes  an  Intermediate  stage  between  a  one-dimensional  auto- regressive  model  described 
earlier  and  a  full  three-dimensional  hybrid  structure  model  to  be  described  In  a  subsequent 
report  The  hybrid  model  has  Its  bases  In  the  one-dimensional  AR  model  and  Is  suggested 
naturally  by  the  present  work. 

2.  THEORETICAL  DISCUSSION 


This  report  is  a  two-dimensional  extension  of  a  previous  report2  that  treated  one-dlmenslonal 
simulation  from  the  autoregressive  (AR)  modeling  perspective.  Though  the  theoretical  background 
for  one-dlmenslonal  analysis  can  be  extended  to  two  dimensions,  the  method  for  estimating  the 
ARMA  coefficients  Is  treated  somewhat  differently  in  this  report. 

For  the  two-dimensional  case,  we  wish  to  simulate  a  horizontal  sheet  of  stationary  data 
having  a  constant  correlation  length  (Lj.),  variance  (o2),  and  spatial  spectra  characterized  by  the 
symmetric  two-dimensional  power  spectral  density  function  (PSD): 


F(fk)  = 


o2va2v 

«(a  a  +  fkT‘ 


’sharma.  R.D..  Duff,  J.W.,  Sundberg.  R.L.,  Gruninger,  J.H..  Bernstein,  L.S.,  Robertson,  D.C.,  and  Healey,  RJ.. 
(1991)  H(gh  Altitude  Atmospheric  Radiance  Code,  Phillips  Laboratory  technical 

2Brown,  J.H..  (1993)  Atmospheric  Structure  Simulation :  An  Autoregressive  Model  for  Smooth  Geophysical  Rower 
Spectra  with  Known  Autocorrelation  Function,  Phillips  Laboratory  technical  report.  PL-TR-93-2185.  ERP  #1128.  ADA27669 1 
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where. 


2V^Hv) 


Lc  and  where,  =  f *  +  f* 


Here,  ffc  is  the  spatial  frequency  (km*1)  and  -2(v+l)  a  the  spectral  slope  of  log(F)  vs  lagfffc). 
The  two-dimensional  ARMA  power  spectral  density  model  is3*4: 

nfr  t 

Parma "  T1Tapw[ 


B(f,.fy) 

A(f„fy) 


where  the  autoregressive  factor  is, 
A(fx'y  )> -  2Ia(m.n),-a— ■T'e-3’"-T> 


m  n 


and  the  moving  average  factor  is, 

m  n 

Tj  and  T2  are  the  sampling  intervals.  Pw  is  the  variance  of  the  white  noise  process,  and  a(m.n), 
b(m.n)  are  respectively  the  autoregressive  and  moving  average  filter  coefficients. 

The  discrete  two-dimensional  series.  x(k.l),  that  approximates  these  deterministic  and  sto¬ 
chastic  processes  can  be  simulated  by  the  finite  filter  difference  equation: 

x(k.l)  =  -  £  ]jT  a(m,n)xfk  -  m.l- nj+ b(m,n)e[k  -  m.l  -  nj 

m  n  m  n 

(nwiMO.O) 


in  which  the  range  of  summations  denote  the  order  of  the  2D  difference  equation.  In  this  report, 
the  order  of  the  arrays  a(m,n)  and  b(m,n)  are  (2M+ 1)  x  (2N+1)  where  the  region  of  support  is  the 
square  full  plane  covering  -M£m£M  and  -N£n£N  and  where  M»N.  The  output  sequence  is  x(k,l), 
and  £(k,I)  is  a  white  noise  input  driving  sequence.  The  a(m,n)  autoregressive  coefficients  were 
calculated  from  an  auto-regressive  estimator  and  the  b(m,n)  moving  average  coefficients  were 
calculated  from  the  equivalent  MA  filter.5 


3Maiple,  S.L.  (1987)  Digital  Spectral  Analysis  with  Applications,  Chapter  6,  Prentice-Hall.  Englewood  Cttfla.  New 
Jersey. 

4 Kay.  Steven  M..(1988)  Modem  Spectral  Estimation,  Theory  A  Application.  Prentice-Hall.  Englewood  Cttfla,  New 
Jersey- 

sUm.  J.S.  (1990)  Two-Dimensional  Signal  and  Image  Processing,  Prentice-Hall,  Englewood  Cttflb  New  Jersey,  pg  269. 


3.  FILTER  CALCULATIONS 


'<et  us  first  treat  the  autoregressive  factor  in  a  region  of  full  plane  support. 


m«-M n»-N 


Considering  the  8-fold  symmetries  of  the  present  problem,  the  following  relationships  apply: 
a(m.n)  *  a(-m,n)  =  a(m.-n)  =  a(-m,-n) 

=  a(n.m)  *  a(-n.m)  =  a(n,-m)  *  a(-n,-m), 
with  similar  symmetries  for  the  b's.  And  we  define  a(0,0)  ■  1. 

Due  to  these  symmetries,  A(fx,fy),  and  in  a  similar  way  B{fx,fy),  can  be  written  as: 


A(f,.f,)=  I;  X  £aio.„* 

m«ln--N  *•  '  n—N 


2«tfrnT, 


e1*  +e-t* 


Using  the  identity:  cos(x)  - - - - ,  then. 


A(fx,fy)  =  2^  ^a(m.n)cos(2jtfxmT1  }e~2*tf,rfTl  +  JalO.nJe-2^^’ 


m=ln=-N 


and  then  similarly. 


A(fx.fy)=4  ^^a(m,n)cos(2xfxnfr1)cos(2xfynT1) 


m>ln«l 


H  « 

+2]£a(0.n)cos(2KfynT, ) + 2^a(m.0)cos(2xfxnfr, )+ a(0.0) 


Assume  that  A(fx,fy)  defines  a  stable  filter,  then  A(fx,fy)  >  0  for  all  (fx.fy);  or. 


/PSD(fx,fy)  =  ^pTjTj 


B(fx,fy) 

A(fx,fy) 


We  then  find  an  approximate  solution  for  A(fx,fy)  by  setting  B(fx.fy)=  . . .  ■.  Then  with 

VPT,T2 

a(0.0)  =  1,  we  find  a(m,n)  for  m£0  and  n£0  by  solving  the  following  least  square  problem  for  a 
large  set  of  fx,fy  values:  1 «  A(fx,fy)^F(fk) .  The  set  of  "theoretical"  PSD  s  used  in  solving  the  least 

square  problem  are  chosen  from  50  circles  defined  by  fk  =  ^fx  +fy  =  constant ,  which  is  equally 

spaced  in  Log^fx  +  fy  and  by  20  points  on  each  circle  defined  by  20  equally  spaced  angles 

tan between  0ant*^.  The  minimum  frequency,  (f,,)^,  used  in  the  fit  is 
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(fk)  =  — — - .  where  N  is  usually  set  to  128.  In  finding  a(m.n)  from  the  least  squares 

nan  N  •  spacing 

procedure,  the  analysis  places  a  weight  of  lOOx  for  the  minimum  frequency. 

Having  found  a  least  squares  solution  for  a(m.n),  we  proceed  as  before  to  solve  for 

M  N 

B(fx.fy)  =  4  ^  b(m,  n)cos(2nfxmTi  )cos(2nfynT1 ) 

n>In>l 

N  M 

+2^bC0.n)cos{2irfynTt)+ 2]^b{m.0)cos(2irfxmTj)+ b(0,0) 

n»l  m=l 


B(f  f  jJpX  T 

by  setting  up  the  least  squares  problem:  1 »  j  1  2  ,  where  A(fx,fy)  are  the  values  found 

from  using  the  solved  set  of  a(m,n)  values,  and  where  the  solution  Is  over  a  large  set  of  fx,fy 
values. 


4.  MODEL  DISCUSSION  AND  RE8ULT8 


The  following  discussion  is  aimed  at  producing  a  practical  two  dimensional  Isotropic  autore¬ 
gressive  /  moving  average  (ARMA)  model  of  horizontal  atmospheric  structure  consistent  with 
having  a  pie-assigned  power  spectral  density.  Correlation  lengths  and  o2  variances  are  taken 
from  Strugala.  et  al.6  In  the  following  discussion  the  full  plane  autoregressive  a(m.n)  coefficients 
will  be  referred  to  as  the  AR  coefficients  and  the  moving  average  b(m.n)  coefficients  will  be  referred 
to  as  the  MA  coefficients.  The  order  of  the  full  plane  AR  or  MA  calculations  will  be  referred  to  as 
having  a  (2M+1)  x  (2M+1)  region  of  support.  Because  the  region  of  support  has  an  8-fold 
symmetry  and  the  desired  2-D  PSD  is  circularly  symmetric,  the  total  number  of  frill  plane  unique 

AR  or  MA  coefficients  is  equal  to  tM  +  2XM+2}i;  whereas  the  total  number  (non-unique)  AR 

coefficients  Is  equal  to  4M(M  + 1)  and  the  total  number  of  (non-unique)  MA  coefficients  is  equal  to 
(2M  +  1)2. 


Three  model  curves  are  plotted  In  each  of  Figures  1-2,  5-6,  9-11,  and  14.  In  these  figures, 
curves  marked  by  an  asterisk  (*)  at  the  first  frequency  are  the  desired  or  "theoretical"  PSD's 
(that  Is,  Ftffc))-  Curves  marked  by  open  squares  (□)  represent  the  PSD  of  the  full  plane  auto- 

(  i  j  la's 


regressive  (AR|  predictor  model  alone 


that  Is, 


A(fx,fy) 


multiplied  by  the  factor  indicated  on 


the  plot  Curves  marked  by  an  open  triangle  (A)  at  the  first  frequency  represent  the  PSD's  of  the 


®Strugala,  LA.,  Newt  J.E..  Futterman.  W..  Schweitzer,  E.L,  Herman,  B.J.,  and  Sean,  R.D.(1991)  Development  of 
High  Resolution  Statistically  Non-stationary  Infrared  Earthltmb  Radiance  Scenes.  Characterization,  Propagation,  and 
Simulation  qf  Sources  and  Backgrounds.  Proceedings  SPtE  -  The  International  Society  Jar  Optical  Engineering.  V1486.  pp 
176-187,  April  1991,  Orlando,  Florida. 
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full  ARMA  predictor  model 


that  ia. 


pT,Ta 


B(f„fy) 


A(fx,fy) 


JJ 


,  which  is  to  say,  the  AR  model  corrected  by 


the  moving  average  |MA]  model).  Except  where  otherwise  noted,  the  sampling  spacing  is  100  m,  N 
for  determination  of  the  minimum  frequency  =  128,  and  the  Nyqulst  frequency  is  5  km'1.  This 
means  that  the  "minimum"  frequency  that  went  into  calculating  the  AR  and  MA  coefficients  was 


f^  =  — •  0.078  km-1.  In  the  model  plots,  this  point  is  marked  by  the  left-most  X.  Figure  1  is 
typical  of  the  model  plots  in  this  report  showing  log- log  power  spectral  densities  (PSD' a).  This  and 


subsequent  two-dimensional  power  spectral  density  plots  have  PSD's  measured  in 


(ST/Temperature)2 

Wavenumber2 


and  wavenumber  measured  in  km'1.  The  input  parameters  for  Figure  1  are 


32  km,  o2  =  0.02,  and  spectral  slope  (S)  «  -8/3.  This  plot  was  calculated  along  a  principal 
axis  using  a  3  x  3  order  AR  calculation  and  a  7  x  7  order  MA  calculation.  Except  for  a  small 
amount  of  overshoot  at  low  frequencies  («  0.2-0. 3  km'1),  agreement  between  desired  and 
predicted  PSD  is  good.  Since  the  ARMA  PSD  model  can  be  calculated  at  any  frequency,  this  and 
subsequent  figures  show  the  expected  bias  between  the  desired  and  predicted  PSD  below  the 
"minimum"  frequency.  Figure  2.  which  is  evaluated  along  a  diagonal  axis,  shows  similar  results 
and  verifies  the  circular  symmetry  of  the  2D  ARMA  PSD.  Further  evidence  of  the  near  circular 
symmetry  is  shown  in  Figure  3.  which  displays  contours  of  constant  PSD  in  two-dimensional 
linear  frequency  space.  Although  perfectly  symmetrical,  the  contours  exhibit  deviations  from 
perfect  circularity.  For  the  most  part  these  deviations  are  suppressed  when  plotted  in  log-log 
frequency  space  as  shown  in  Figure  4. 


Figures  5  and  6  use  the  same  parameters  as  Figure  2,  except  that  Figure  5  is  evaluated  for  a 
slope  of  -3  and  Figure  6  is  evaluated  for  a  slope  of  -11/3.  These  plots  indicate  that  slightly  better 
agreement  obtains  between  the  "theoretical”  and  model  curves  as  the  spectral  slope  steepens.  Also 
the  "DC”  bias  diminishes  with  steeper  slope.  Figures  7  and  8  show  constant  contours  of  the  -11/3 
model  PSD  plotted  in  two-dimensional  linear  space  and  two-dimensional  log-log  space 
respectively.  Comparing  these  to  Figures  3  and  4,  we  see  that  slightly  better  circularity  evidently 
obtains  at  the  steeper  slope. 


Figures  9,  10.  and  1 1  repeat  the  input  parameters  of  Figures  2,  5,  and  6.  except  that  is 
increased  to  84  km  and  o2  is  decreased  to  0.005.  Again  slightly  better  agreement  obtains  as  the 
slope  steepens  from  -8/3  in  Figure  9  to  -3  in  Figure  10  to  -1 1/3  in  Figure  1 1.  In  fact,  agreement 
is  particularly  good  at  the  "minimum"  frequency  for  the  -11/3  slope.  Figures  12  and  13  which 
show,  respectively,  two-dimensional  linear  and  log-log  frequency  plots  of  constant  PSD  contours 
for  =  84  km  and  -11/3  spectral  slope,  are  comparable  to  the  =  32  km  contour  plots  of 
Figures  7  and  8. 


Figure  14  repeats  the  calculation  of  Figure  1  for  =  32  km,  2-D  slope  =  -8/3,  but  now  N  * 
1024.  Clearly  the  low  order  ARMA  model  does  not  provide  good  agreement  for  the  wider  frequency 
range.  Figure  15  provides  an  indication  of  the  results  of  using  higher  order  AR  and  MA  models. 
The  curves  in  Figure  15  also  are  calculated  for  =  32  km,  2-D  slope  =  -8/3,  and  N  *  1024,  and 
thus  may  be  compared  with  Figure  14.  Comparing  Just  the  AR  curves,  we  observe  that  as  the 
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order  of  the  AR  process  increases  from  3x3  (Figure  14)  to  13  x  13  (Figure  15):  the  low  frequency 
"rollover"  occurs  at  smaller  frequencies,  and  the  number  of  "ripples"  Increase.  Comparing  the 
complete  ARMA  curves,  we  observe  that  as  the  order  of  the  AR  process  increases  from  3  x  3  to  13 
x  13  and  as  the  order  of  the  MA  process  increases  from  7  x  7  to  25  x  25:  the  "theoretical”  and 
model  PSD's  agree  much  better  at  low  frequency,  the  "DC"  bias  decreases,  and  the  amplitude  of 
the  "ripples”  decrease.  As  observed,  PSD's  having  a  small  slope,  large  correlation  length,  and 
small  data  spacing  are  more  difficult  to  model  than  PSD's  with  larger  slope,  smaller  correlation 
length,  and  wider  data  spacing. 

The  analysis  indicates  that  two-dimensional  Isotropic  horizontal  atmospheric  temperature 
structure,  having  a  data  spacing  of  100  m  and  a  characteristically  smooth  power  spectral  density 
can  be  modeled  within  a  frequency  band  ranging  from  0.08  km*1  to  5  km*1  by  an  autoregressive  3 
x  3  order  process  and  a  moving  average  7x7  order  process.  At  the  expense  of  increasing 
computer  time,  one  may  increase  the  frequency  range  and  achieve  agreement  to  0.01  km*1  by 
choosing  a  13  x  13  AR  model  and  25  x  25  MA  model. 


7 


*•  Ull-B  MUM  LM»  na  art  U.4* 

*  Nn  *  «  S  Ml  NU  7  17 

■in.  nv  .i/iiMki.  7»U-M  k#cnc  H.ui  ms  a*  i.ok-oi 


10® 


104 


IO3 


10* 


101 
l  io° 

s 

V) 

i  10* 1 

-I 

t 

|  .0-2 
10-3 


10-4 


io-s 


io-s 


io-7 

10'S  IO-4  10"3  10'2  IO-1  10°  10* 

MPVENUMERI  1./MH1 

Figure  1.  Log-log  {dot  of  two-dimensional  power  spectral  density  (PSD)  versus  spatial  frequency 
{dotted  along  a  major  (North -South  or  East-West)  axis.  Curve  marked  by  an  asterisk  (•) 
at  the  first  frequency  Is  the  desired  or  "theoretical"  PSD.  Curve  marked  by  open  squares 
(o)  represents  the  PSD  of  the  full  plane  auto-regressive  |AR]  predictor  model  multiplied 
by  the  factor  indicated  on  the  plot  Curve  marked  by  an  open  triangle  (A)  at  the  first 
frequency  represents  the  PSD  of  the  full  ARMA  predictor  model  (which  Is  to  say.  the  AR 
model  corrected  by  the  moving  average  [MAI  model;  this  model  is  found  by  weighted 

least  squares  fit  where  the  minimum  frequency  of  the  fit  is  - - - - . ).  The  2-D  slope 

N*  spacing 

is  -2.67  and  the  correlation  length  >  32  km.  The  spacing  in  this  and  all  subsequent 
PSD'a  *  100  m.  Order  of  AR  coefficients  >3x3.  Order  of  MA  coefficients  >7x7. 
N  >  128  in  this  and  other  plots  unless  otherwise  noted 
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Figure  2.  Same  as  Figure  1,  but  the  spatial  frequency  is  (dotted  along  a  diagonal  axis 
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Figure  4.  Same  as  Figure  3  but  the  two-dimensional  plot  Is  (dotted  in  logarithmic  frequency  space 
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Figure  7.  Same  as  Figure  3;  1^-32  km.  except  2-D  slope  ■  -3.67 
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Figure  8.  Same  as  Figure  4;  L*.  =  32  km.  except  2-D  slope  =  -3.67 
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Figure  9.  Same  as  Figure  2.  except  Lc  «  84  km  and  2-D  slope  « 


-2.67 
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Figure  10.  Same  as  Figure  2,  except  1^-84  km  and  2-D  slope  ■  -3 
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Figure  1 1.  Same  as  Figure  2.  except  Lg  «  84  km  and  2-D  slope  ■  -3.67 
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Figure  12.  Same  as  Figure  3,  except  Lg  ■  84  km  and  2-D  slope  ■  -3.67 
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Figure  13.  Same  as  Figure  4.  except  ■  84  km  and  2-D  slope  ■  -3.67 
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Figure  14.  Same  as  Figure  1.  except  number  of  frequencies.  N  »  1024 
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Figure  15.  Same  as  Figure  14.  except  the  order  of  AR  coefficients  -  13  x  13.  Older  of  MA 
coefficients  «  25  x  25,  N  ■  1024. 


22 


8.  SIMULATION  •  PROCEDURE 


To  simulate  sets  of  two-dimensional  Gaussian  correlated  data  from  the  model  filters  described 
above,  we  start  out  with  a  sheet  of  pseudo- random  numbers  having  a  standard  deviation  o  »  1. 
The  square  N  x  M  sheet  (N  ■  M  »  1300)  Is  somewhat  larger  than  the  sheet  we  wish  to  simulate, 
say  an  extra  138  values  on  each  side.  The  extra  points  allow  iterative  relaxation.  Figure  16  depicts 
the  1300  x  1300  two-dimensional  sheet  of  random  numbers. 

A  problem  ensues  when  we  attempt  to  employ  the  full  plane  A(fx,fy)  filter.  A(fx,fy)  is  not  a 
"computable"  Infinite  impulse  response  filter.  A  "computable"  filter  can  be  applied  directly  to  a 
sheet  of  numbers.  Affx,fy),  however,  involves  the  solution  of  N2  linear  equations  In  N2  unknowns. 
For  N  -  1000,  the  computational  requirements  become  prohibitive.  An  approximate  solution  Is 
therefore  obtained  as  follows. 

An  Initial  "guess”  for  the  solution  of  the  N2  equations  In  N2  unknowns  Is  made  by  summing 
the  separate  results  of  four  quarter  plane  filters,  AAi(fx,fy).  applied  to  each  comer  of  the  same 
sheet  of  random  numbers.  Each  of  the  quarter  plane  filters  is  "computable".  In  effect,  computable 
means  that  the  resulting  N2  linear  equations  in  N2  unknowns  can  be  written  as  the  solution  of  a 
triangular  matrix  equation.  The  four  quarter  plane  filters  are  derived  from  the  full  plane  filter 
coefficients  as  follows,  where  a(m,n)  are  the  full  plane  AR  filter  coefficients  and  al(m.n),  a2(m,n), 
a3(m,n),  and  a4(m.n)  are  the  four  quarter  plane  filter  coefficients: 


Let 

al(m,n)  =  0 

for 

m<0or 

n<0 

al(m,n)  =  2a(m,n) 

m  =  0. 

n>0 

al(m,n)  =  2a(m,n) 

m  >0, 

n  =  0 

al(m,n)  =  4a(m,n) 

m>0. 

n  >0 

Let 

a2(m,n)  =  0 

for 

m>0or 

n<0 

a2(m.n)  =  2a(m,n) 

m  =  0, 

n>0 

a2(m.n)  =  2a(m,n) 

m<0. 

n  =  0 

a2(m,n)  =  4a(m,n) 

m<0. 

n  >  0 

Let 

a3(m,n)  =  0 

for 

m  <0or 

n>0 

a3(m,n)  =  2a(m,n) 

m  =  0, 

n  <  0 

a3(m,n)  =  2a(m,n) 

m  >0, 

n  =  0 

a3(m,n)  =  4a(m,n) 

m  >0, 

n  <0 

Let 

a4(m,n)  =  0 

for 

m  >0or 

n>0 

a4(m,n)  =  2a(m,n) 

m  =  0. 

n<0 

a4(m,n)  =  2a(m,n) 

m  <0, 

n  =  0 

a4(m,n)  =  4a(m,n) 

m<0. 

n  <0 

23 


24 


Figure  16.  Pictorial  representation  of  two-dimensional  sheet  of  Gaussian  random  numbers 


These  relationships  are  schematically  Illustrated  In  Figure  17.  The  principal  reason  for  employing 
this  Initialization  procedure  is  that  It  seems  to  allow  the  iteration  scheme  described  below  to 
converge  faster  to  the  actual  solution.  In  general,  Iteration  schemes  that  attempt  to  simulate 
filtered  data  converge  Inversely  with  the  power  spectral  density  as  a  function  of  frequency.  Since 
the  greatest  power  occurs  at  the  lowest  frequencies,  we  therefore  wish  to  devise  an  initialization 
procedure  that  enables  the  low  frequency  behavior  of  the  simulation  to  follow  the  low  frequency 
behavior  of  the  AR  filter.  This  is  accomplished  with  near  circular  symmetry  by  summing  the 
results  of  the  four  quarter  plane  filters.  Figures  18,  19.  and  20  Illustrate  the  method  for  applying 
the  quarter  plane  filter.  Figure  18  depicts  the  starting  condition  for  a  2  x  2  quarter  plane  AR 
example  (derived  as  Implied  above  from  a  3  x  3  full  plane  AR  model;  that  is,  O  q.p. 
«  (O  f.p.  +  1}  /2).  The  procedure  Is  iterative  and  starts  In  the  upper  right  comer.  The  arrows 
indicate  that  the  filter  is  moved  column-wise  to  the  left.  At  the  left  edge  of  the  sheet,  the  filter  Is 
brought  back  to  the  right  and  dropped  one  row.  Figure  19  shows  an  intermediate  stage  in  the  first 
comer  process.  It  illustrates  that  points  within  and  to  the  right  of  the  filter  have  been  processed, 
while  points  in  the  adjacent  first  row  to  the  left  of  the  filter  region  also  have  been  processed.  In 
this  particular  iteration,  the  point  at  the  extreme  left-bottom  of  the  2x2  filter  region  contains  the 
original  random  number  but  gets  replaced  by  the  2x2  AR  filtered  value.  Note  that  points  in  the 
one-point-wide  border  remain  unprocessed.  Figure  20  Illustrates  the  starting  condition  for  the 
second  pass  at  the  upper  left  comer.  Here  the  filter  is  moved  to  the  right  and  then  down.  The 
extreme  right-bottom  point  In  the  filter  region  contains  the  original  random  number  but  gets 
replaced  by  the  2  x  2  filtered  value.  After  four  such  passes,  (one  pass  starting  from  each  comer 
and  each  pass  using  the  same  original  random  numbers)  the  results  are  summed.  The  following 
analysis  indicates  why  the  four-quarter-plane  solution  provides  a  good  approximation  of  the  low 
frequency  behavior. 

We  have. 

1  1  1  1  1 

AA1(fx,fy)  AA2(fx.fy)  AA3(fx,fy)  AA4(fx.fy) 

where  AA(fx,fy)  is  the  AR  filter  associated  with  the  sum  of  the  four  quarter  plane  filters  and 
AAj(fx.fy)  is  the  AR  filter  associated  with  the  i^1  quarter  plane  filter. 

Using  the  identity,  e-lx  =  cos(x)  -  isin(x)  and  setting  N  *  M  and  Tj  ■  T2,  we  have, 

1  _ 1 _ 

AAi(fx-fy)  ^  ^a,(m.n){cos[2xT,(mfx  +nfy)]-lsin[2nT1(mfx  +nfy)]} 

m«-M  n*-N 

and  as  fx  -» 0  and  fy  -» 0 

_ 1 _  _ 1 _ 

AAi(fx  ->0.fy  -»0)  ^  ^a,(m.n)cos(2xT1mfx)cos(2xT,nfy) 
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Figure  18.  Diagram  of  starting  condition  for  application  of  first  quarter  plane  AR  Alter.  Example  depicted  for  2  x  2  AR  quarter 
plane  region  of  support 
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Figure  19.  Diagram  of  Intermediate  step  In  processing  first  of  four  quarter  plane  In  itlallzatkms.  2  x  2  AR  quarter  plane  example 


Ion  for  application  of  second  quarter  plane  AR  filter.  Depicted  for  same  2  x  2  AR  quarter 


1 


and.  noting  that  aj(m,n)  ■  a2(-m,n)  *  agt-m.-n)  *  a^m.-n), 


m  n  ■  n 

^  ^al(m,n)co8(2*T1mfx)co8(2jtT1nfy)  ^^aI(m.n)co6(-2xT1mf1)cos(-2icT1nfy) 


m-0n-0 


m-On-0 


£  ^Tal(m,  n)cos(-2*T1mfx  )cos(2xT,nfy )  £  £  al(m,  n)cos(2<r1mfz  )cos(-2ifl,,nfy ) 

m*0n«0  m«0n«0 


m  n 

££al(m.n)cos[2*T,(mfx  +  nfy)] 


m*0n~0 


m  m 

£  £  al(m.  n)[cos(2ionfxT,  )cos(2jmfyT, )  -  sin(2xmfxT1  )sin(2xnfyT, )] 


m-On-O 


m  n 

£  £al(m.  n)[cos(2xinfxT,  )cos(2xnfyT1 )] 


m-On-O 


where  again,  we  have  Ignored  the  sine  terms.  Recalling  the  symmetries  of  the  aj(m.n)  coefficients, 
we  find, 

M  N 

4AA(fx.fy)  =  4  y^a(m,  n)cos(2n  fxmT,)cos(2x  fynT,) 

m°ln*l 

N  M 

+  2Ta(0,n)cos(2itfynTl)  +  2  ^a(m.0)cos(2n  fxmT,)  +  a(0,0) 

n»l  m>l 

1  4 

and  arrive  at  the  approximation.  — —  -  —  It  Is  unclear  why  the  sum  of  the  filters 

AaUx.V  A(fx,fy) 

fj-tO 

fy-»0 

performs  better  than  the  average.  A  more  complete  description  of  this  process  for  the  four  quarter 
plane  filters  is  provided  in  the  appendix. 

The  four-quarter-plane  procedure  provides  an  Initial  guess  for  the  correct  data  array. 
Subsequently,  this  initial  guess  Is  used  as  input  to  the  full  plane  process.  Following  the 
four-quarter-plane  Initialization,  the  procedure  continues  by  iteratively  updating  and  correcting 
the  spatial  data  set  In  the  following  way.  At  a  point  in  the  array  (x,y)  we  sum  the  products  of  all 
neighboring  points  within  the  region  of  support  and  their  corresponding  full  plane  coefficients.  Tb 
this  sum  Is  added  the  original  (x.y)  random  number  e(x.y);  the  result  replacing  the  latest  (x.y) 
value.  This  procedure  continues  for  each  point  (x.y)  in  the  array.  Array  substitutions  are  made  In 
the  same  order  as  taken  for  the  four  quarter  plane  filters.  To  vector  optimize  for  the  Convex 
computer,  elements  in  the  lower  row  of  the  filter  region  were  not  replaced  until  all  columns  of  that 
row  were  processed.  The  optimization  process  decreased  computer  run  times  significantly.  One 
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set  of  substitutions  (a  set  being  defined  by  the  same  and  every  path  taken  by  the  four  passes  as 
previously  defined  for  the  four  quarter  plane  passes  -  one  pass  per  comer)  constitutes  one  total 
Iteration  In  this  stage  of  the  procedure. 

Figures  21  and  22  Illustrate  application  of  the  full  plane  autoregressive  filter  A(fx,fy)  using 
multiple  predictor-corrector  iterations.  Figure  21  shows  the  starting  condition  for  a  3  x  3  full 
plane  AR  example.  The  plane  is  initialized  by  the  final  (summed)  values  of  the  quarter  plane  AR 
process  as  described  above.  As  before,  the  procedure  is  iterative  and  starts  in  the  upper  right 
comer.  Also,  as  described  before,  points  in  a  one-point-wide  border  remain  unprocessed.  Figure 
22  shows  an  intermediate  stage  in  the  first  comer  process.  It  illustrates  that  points  above  the 
center  point  and  points  to  the  right  in  the  row  adjacent  to  the  center  point  of  the  filter  have  been 
processed.  In  this  iteration,  the  center  point  of  the  3x3  filter  region  contains  the  original  random 
number  but  gets  replaced  by  the  3x3  full  plane  AR  filtered  value.  Pour  such  passes,  one  starting 
from  each  comer,  are  performed.  Unlike  the  quarter  plane  initialization  process,  each  subsequent 
pass  is  performed  on  the  numerical  results  of  the  preceding  pass.  Consequently,  no  summing  is 
performed.  We  emphasize,  however,  that  every  iteration  begins  with  a  center  value  equal  to  the 
original  random  number,  which  gets  replaced  by  the  full  plane  AR  filtered  value.  As  discussed 
earlier,  it  Is  necessary  to  cycle  though  the  full  plane  AR  procedure  several  times  to  achieve  a 
satisfactory  estimate  of  the  set  of  1300  linear  equations  In  1300  unknowns.  In  general,  we  find 
that  20  cycles  provides  reasonable  power  spectral  densities  that  parallel  the  "theoretical''  PSD's. 

Figure  23  Illustrates  application  of  the  final  moving  average  filter  B(fx,fy).  The  MA  process  Is  a 
mapping  procedure  that  Is  performed  only  once  on  the  sheet  containing  the  final  numerical 
results  of  the  full  plane  process.  As  before,  we  start  from  the  upper  right  comer  and  proceed  left 
and  down.  Actually  we  could  start  and  end  anywhere  since  the  process  simply  maps  the  center 
points  of  the  7x7  MA  full  plane  region  to  points  co-located  on  a  separate  sheet 

6.  ARMA  SIMULATIONS  -  RESULTS  AND  DISCUSSION 

In  the  above  procedure  we  discussed  a  method  for  generating  simulated  two-dimensional 
correlated  data  by  applying  two-dimensional  ARMA  filters  to  Gaussian  random  noise.  We  now 
examine  several  simulated  data  sets  for  fidelity  to  the  original  or  "theoretical"  power  spectral 
density  specifications. 

Figure  24  shows  four  logarithm!  plots  of  two-dimensional  power  spectral  density  (PSD) 
versus  logarithmic  spatial  frequency  plotted  along  a  diagonal  axis.  As  before,  the  curve  marked  by 
an  asterisk  (*)  Is  the  desired  or  "theoretical"  PSD.  The  curve  marked  by  open  squares  (□) 

represents  the  PSD  of  the  full  plane  auto-regressive  (AR]  predictor  model.  The  curve  marked  by  an 
open  triangle  (A)  at  the  first  frequency  represents  the  PSD  of  the  full  ARMA  predictor  model  (which 
is  to  say.  the  (AR]  model  corrected  by  the  moving  average  [MA]  model).  Finally,  the  curve  marked 
by  an  X  is  a  two-dimensional  perlodogram  derived  from  two-dimensional  simulated  data.  The 
1024  x  1024  simulated  data  set  was  derived  from  an  original  set  of  1300  x  1300  random 
numbers.  Figure  24  is  calculated  for  a  "theoretical''  2-D  slope  of  -8/3  and  correlation  length  = 
32  km.  The  region  of  support  for  the  full  plane  AR  model  is  of  order  3x3  and  the  region  of 
support  for  the  full  plane  MA  model  is  of  order  7x7. 
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Figure  21.  Diagram  of  stating  condition  for  application  of  first  full  plane  AR  filter.  Example  depicted  for  3  x  3  AR  full  plane  region  of  support 
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Figure  22.  Diagram  of  intermediate  step  in  processing  first  of  4  x  I  foil  plane  filters  (I  *  20  is  thr  number  of  times 
Depicted  for  3  x  3  AR  full  plane  region  of  support.  Filter  uses  original  random  number  for  center  point 
surrounding  points.  Center  point  gets  replaced  by  filter  value 
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Figure  24.  Log-log  plot  of  two-dimensional  power  spectral  density  (PSD)  versus  spatial  frequency 
plotted  along  a  diagonal.  Curve  marked  by  an  asterisk  (•)  is  the  desired  or  "theoretical" 
PSD.  Curve  marked  by  open  squares  (□)  represents  the  PSD  of  the  full  plane 
auto-regressive  [AR]  predictor  model  multiplied  by  the  factor  indicated  on  the  {dot 
Curve  marked  by  an  open  triangle  (A)  at  the  first  frequency  represents  the  PSD  of  the 
full  ARMA  predictor  model  (which  Is  to  say,  the  AR  model  corrected  by  the  moving 
average  [MAI  model).  Curve  marked  by  an  X  is  the  2-D  pertodogram  derived  from  the 
2-D  simulation.  The  1024  x  1024  simulated  data  set  was  derived  from  an  original  set 
of  1300  x  1300  random  numbers  (see  text  for  discussion  of  DFT  low  frequency 
replacement).  The  "theoretical"  2-D  slope  is  -2.67  and  the  correlation  length  * 
32  km.  The  order  of  the  AR  coefficients  >3x3.  Order  of  the  MA  coefficients  >7x7 
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As  noted,  the  2-D  FFT  periodogram  corresponding  to  Figure  24  was  calculated  from  the  1024 
x  1024  simulated  data  array.  However,  since  a  goal  of  this  effort  was  to  mintmtee  the 
computational  burden,  the  minimum  frequency  that  was  adequately  modeled  by  the  3x3  AR  and 

7x7  MA  processes  occurred  at  a  minimum  frequency  of  L*.  =  — —  - - - =  0.078  km-1.  To 

NAx  128x0.1 

expand  the  usable  frequency  range  and  lower  the  minimum  simulated  frequency  to 

fman  =  —  -  1  =  0.0098  km"1  the  following  procedure  was  performed.  At  each  frequency  below 

1024x0.1 

0.078  km'1,  the  power  (as  determined  by  the  2-D  FFT  periodogram)  was  subtracted  from  the  1024 
x  1024  simulated  data  array.  A  short  12  x  12  discrete  Fourier  transform,  DFT,  procedure  was 
then  performed  to  replace  the  simulated  power  from  0.0098  km*1  to  0.078  km1.  A  detailed 
description  of  the  procedure  may  be  found  In  the  appendix.  Agreement  between  the  2-D  FFT 
periodogram  and  the  ARMA  model  Is  good  to  0.078  km*1.  With  the  addition  of  the  DFT 
replacement  procedure,  agreement  with  the  "theoretical”  PSD  curve  is  good  over  the  entire 
frequency  range.  0.0098  km*1  to  5  km*1. 

One  may  also  examine  the  simulation  by  comparing  the  "theoretical”  one-dimensional  power 
spectral  density  of  the  two-dimensional  data  to  one-dimensional  equivalents  of  the  simulated 
data.  Figure  25  shows  three  logarithmic  plots  of  one-dimensional  power  spectral  density  (PSD) 
versus  logarithmic  spatial  frequency  plotted  along  a  major  axis.  The  curve  marked  by  an  asterisk 
(*)  is  the  desired  or  "theoretical”  1-D  PSD.  The  curve  marked  by  an  open  circle  (o)  at  the  first 
frequency  represents  the  1-D  PSD  derived  by  integrating  the  2-D  PSD  ARMA  model.  Finally,  the 
curve  marked  by  an  X  at  the  first  frequency  represents  a  one-dimensional  autoregressive  PSD 
model  derived  from  the  simulated  2-D  data.  The  latter  is  accomplished  by  treating  all  rows  and 
columns  as  1-D  data  and  finding  the  corresponding  AR  model.  The  parameters  are  the  same  as 
used  In  Figure  24  with  the  2-D  -8/3  slope  becoming  a  1-D  slope  of  -5/3.  Figure  25  shows  that  the 
1-D  power  spectral  density  curve  as  derived  by  Integrating  the  2-D  ARMA  model  is  faithful  to  the 
"theoretical"  model.  Of  course  the  curve  rolls  off  as  it  approaches  its  lower  limit  Also,  the  1-D  AR 
power  spectral  density  curve  follows  the  "theoretical"  to  much  lower  frequencies.  The  slight  over¬ 
shoot  of  the  AR-PSD  model  seems  to  be  due  to  the  difficulty  In  approximating  the  PSD  at  low 
frequencies  using  an  AR  model  of  limited  extent. 

We  may  also  examine  the  one-dimensional  FFT  periodogram  derived  from  the  simulated  two- 
dimensional  data  set  Figure  26  shows  the  same  curves  as  Figure  25,  except  that  the  curve 
marked  by  an  X  Is  the  one-dimensional  power  spectral  density  derived  by  finding  the  FFT  of  all 
rows  and  columns.  Again,  good  agreement  obtains  between  the  "theoretical"  and  the  1-D 
periodogram.  Of  course,  at  the  expense  of  computer  time,  better  and  better  agreement  would 
obtain  (especially  near  the  "rollover"  frequency)  as  the  number  of  AR  iterations  increase. 

Figures  27  and  28  show  the  same  four  two-dimensional  PSD  plots  as  Figure  24  accept  that 
Figure  27  is  plotted  for  a  slope  of  -3  and  Figure  28  is  plotted  for  a  slope  of  -11/3.  These  cases 
show  that  the  simulated  data  retains  high  fidelity  to  the  "theoretical”  2-D  PSD  slope  from  -8/3  to 
-11/3. 
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Figure  25.  Log-log  plot  of  one-dimensional  power  spectral  density  (PSD)  versus  spatial  fin 
plotted  along  a  major  axis.  Curve  marked  by  an  asterisk  (•)  is  the  do 
"theoretical"  1-D  PSD.  Curve  marked  by  an  open  circle  (o)  at  the  first  fin 
represents  the  1-D  PSD  derived  by  Integrating  the  2-D  PSD  ARMA  model 
marked  by  an  X  at  the  first  frequency  represents  an  autoregressive  PSD  model 
from  the  simulated  2-D  data  (with  DFT  low  frequency  replacement).  The  1024 
simulated  data  set  was  derived  from  an  original  set  of  1300  x  1300  random  m 
The  theoretical"  1-D  slope  Is  -1.67  and  the  correlation  length  Lg  ■  32  km.  The 
the  AR  coefficients  >3x3.  Order  of  the  MA  coefficients  >7x7 
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Figure  26.  Log-log  plot  of  one-dimensional  power  spectral  density  (PSD)  versus  spatial  frequency 
(dotted  along  a  major  axis.  Curve  marked  by  an  asterisk  (•)  is  the  desired  or 
"theoretical"  1-D  PSD.  Curve  marked  by  an  open  circle  (o)  at  the  first  frequency 
represents  the  1-D  PSD  derived  by  Integrating  the  2-D  PSD  ARMA  model.  Curve 
marked  by  Xs  represents  the  ID  periodogram  derived  from  the  rows  and  of 

the  2-D  simulated  data  (with  DFT  low  frequency  replacement).  The  1024  x  1024 
simulated  data  set  was  derived  from  an  original  set  of  1300  x  1300  random  numbers. 
The  "theoretical"  1-D  slope  is  -1.67  and  the  correlation  length  >  32  km.  The  order  of 
the  AR  coefficients  >3x3.  Order  of  the  MA  coefficients  >7x7 
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Figure  27.  Same  as  Figure  24,  except  2-D  slope  ■  -3.00 
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Same  as  Figure  24,  except  2-D  slope  ■  -3.67 


The  curves  presented  In  Figures  29.  30.  and  31  repeat  the  plots  of  Figures  24,  27.  ami  28 
except  that  L*.  ■  84  km.  Clearly  the  two-dimensional  periodograms  indicate  good  agreement  with 
the  "theoretical'’  PSD  slopes  for  the  entire  spread  of  correlation  lengths  and  spectral  slopes. 

The  2-D  curves  presented  in  Figure  32  repeat  the  [dot  of  Figure  29  except  Figure  32  is 
calculated  for  an  AR  order  13  x  13,  MA  order  25  x  25,  N  ■  1024,  and  with  no  low  frequency 
replacement.  Clearly,  the  AR  and  ARMA  model  curves  match  the  "theoretical”  PSD  better 
especially,  the  high  order  ARMA  model  matches  the  "theoretical”  PSD  to  much  lower  frequencies. 
Also,  the  perlodogram  matches  the  "theoretical"  PSD  well.  Not  much  improvement,  however,  was 
seen  as  the  AR  order  increased  above  order  13  x  13.  Again,  the  cost  of  Increasing  the  ARMA  model 
to  high  order  is  greatly  increased  computer  time. 

The  2-D  curves  presented  in  Figure  33  repeat  the  plots  of  Figure  29,  except  that  the 
perlodogram  plot  in  this  figure  shows  the  effect  of  subtracting  (Hit  the  low  frequency  components 
from  the  1024  x  1024  ARMA  simulated  scene.  The  PSD  falls  abruptly  by  many  orders  of 
magnitude  just  below  the  128  point  "minimum"  cutoff  frequency.  Figure  34  presents  a  two- 
dimensional  gray  scale  102.4  km  x  102.4  km  pictorial  realization  corresponding  to  Figure  33.  The 
high-frequency-only  structured  data  was  constructed  for  L^,  *  84  km,  2-D  slope  =  -8/3,  spacing  = 
100  m,  AR  order  =  3x3,  and  MA  order  =  7x7  but  with  the  low  frequency  components  removed. 
The  two  figures  show  that  the  procedure  described  in  the  appendix  is  highly  effective  in  removing 
the  low  frequency  components,  leaving  only  the  ARMA  high  frequency  components.  Figure  35 
repeats  the  curves  presented  in  Figure  29.  except  that  the  perlodogram  shows  only  the  added  low 
frequency  DFT  replacement  portion  of  the  spectrum.  Clearly,  the  12  x  12  replacement  scheme, 
described  in  the  appendix,  is  able  to  provide  low  frequency  power  components  that  match  the 
"theoretical"  PSD  well. 

For  completeness.  Figure  36  shows  a  finished  two-dimensional  gray  scale  pictorial  realization 
of  102.4  x  102.4  km  ARMA  simulated  data  with  low  frequency  DFT  replacement.  The  structured 
data  were  constructed  for  L^.  =  84  km.  2-D  slope  *  -8/3.  spacing  =  100  m,  AR  order  =  3x3.  and 
MA  order  =  7x7.  The  model  and  calculated  power  spectral  densities  of  this  data  were  presented 
in  Figure  29. 


41 


I0'5  to"4  JO-3  10_Z  10_1  10°  10* 

UPVENUneERU./Kn)  SaRTIUX«Z+wr*«2) 


Figure  29.  Log-log  plot  of  two-dimensional  power  spectral  density  (PSD)  versus  spatial  frequency 
plotted  along  a  diagonal.  Curve  marked  by  an  asterisk  (•)  is  the  desired  or  "theoretical'’ 
PSD.  Curve  marked  by  open  squares  (o)  represents  the  PSD  of  the  full  plane 
auto- regressive  [AR]  predictor  model  multiplied  by  the  factor  indicated  on  the  plot 
Curve  marked  by  an  open  triangle  (A)  at  the  first  frequency  represents  the  PSD  of  the 
full  ARMA  predictor  model  (which  is  to  say,  the  AR  model  corrected  tty  the  moving 
average  IMA]  model).  Curve  marked  by  an  X  Is  the  2-D  periodogram  derived  from  the 
2-D  simulation.  The  1024  x  1024  simulated  data  set  was  derived  from  an  original  set 
of  1300  x  1300  random  numbers  (see  text  for  discussion  of  DFT  low  frequency 
replacement).  The  "theoretical”  2-D  slope  Is  -2.67  and  the  correlation  length  1^,  * 
84  km.  The  order  of  the  AR  coefficients  >3x3.  Order  of  the  MA  coefficients  *7x7 
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Figure  30.  Same  as  Figure  29,  except  2-D  slope  *  -3 
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Figure  31.  Same  as  Figure  29,  except  2-D  slope 
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Figure  32.  Same  as  Figure  29,  except  the  order  of  the  AR  filter  =  13  x  13  and  the  order  of  the  MA 
filter  =  25x25  (no  DFT  low  frequency  replacement  values) 
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Figure  33.  Same  as  Figure  29,  except  the  low  frequency  power  has  been  subtracted  from  the 
1024  x  1024  simulated  data 
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Figure  34.  Two-dimensional  gray  scale  pictorial  representation  of  1024  x  1024  ARMA  simulated 
data  with  low  frequencies  subtracted  (as  in  Figure  33).  The  data  set  was  constructed 
for  =  84  km.  2-D  slope  =  -8/3.  spacing  *=  100  m,  AR  order  *3x3,  and  MA  order 
■  7x7 
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Figure  35.  Same  as  Figure  29,  except  only  the  DFT  low  frequency  replacement  power  values  are 
shown  In  the  simulation 


Figure  36.  Two-dimensional  gray  scale  pictorial  representation  of  1024  x  1024  ARMA  simulated 
data  (corresponds  to  Figure  29).  The  data  set  was  constructed  for  =  84  km.  2-D 
slope  *  -8/3,  spacing  =  100  m,  AR  order  =  3x3,  and  MA  order  =7x7 
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7.  CONVENTIONAL  SIMULATION  -  RESULTS  AND  DISCUSSION 


This  section  focuses  on  comparing  the  2-D  ARMA  results  with  a  corresponding  aimnb»Hnn 
obtained  by  invoking  the  conventional  fast  Fourier  transform  method.  The  FFT  method  is 
straightforward  in  that  (a)  we  calculate  the  FFT  of  e(x,y),  where  e(x.y)  is  chosen  to  be  the  same  set 
of  Gaussian  random  numbers  that  we  used  in  calculating  the  ARMA  simulation;  (b)  multiply  by 
the  desired  filter  (that  is.  the  square  root  of  the  desired  PSD.  F(fx,f^).  divided  by  the  spacing,  Tj); 
and  (c).  perform  the  inverse  FFT.  Where  g(x.y)  comprises  the  conventional  2-D  simulation,  the 
process  is  described  by: 


g(x.y)  *  lift 


{fft(eCx.yl)} 


x 


T, 


Figure  37  shows  a  two  dimensional  periodogram  of  data  simulated  by  using  this  technique. 
The  log-log  plot  shows  the  two-dimensional  power  spectral  density  (PSD)  versus  spatial  frequency 
plotted  along  a  diagonal.  The  curve  marked  by  an  asterisk  (•)  is  the  desired  or  "theoretical''  PSD 
and  the  curve  marked  by  an  X  is  the  2-D  PSD  obtained  from  the  1024  x  1024  FFT  simulated  data 
set  The  "theoretical''  2-D  slope  is  -2.67  and  the  correlation  length  L*.  *  84  km.  Figure  37  thus 
corresponds  to  the  ARMA  plot  in  Figure  29.  Examination  of  these  plots  show  that  the  ARMA  and 
FFT  simulation  methods  appear  to  yield  very  similar  PSD's. 

Figure  38  shows  a  one-dimensional  periodogram  of  the  two-dimensional  data  set  Just 
described.  The  log-log  plot  shows  the  1-D  power  spectral  density  (PSD)  versus  spatial  frequency 
plotted  along  a  major  axis.  Again,  the  curve  marked  by  an  asterisk  (•)  is  the  desired  or 
"theoretical"  1-D  PSD  and  the  curve  marked  by  Xs  represents  the  ID  periodogram  derived  from 
the  rows  and  columns  of  the  2-D  simulated  data  corresponding  to  Figure  37.  Not  unexpectedly, 
good  agreement  obtains  between  the  "theoretical"  1-D  PSD  and  the  PSD  of  the  FFT  generated  data 
set 

For  completeness.  Figure  39  shows  a  two-dimensional  gray  scale  pictorial  realization  of 
102.4  x  102.4  km  FFT  simulated  data  corresponding  to  Figure  37.  The  structured  data  was 
constructed  for  1^,  =  84  km.  2-D  slope  =  -8/3,  and  spacing  =  100  m.  Comparison  of  the  2-D  ARMA 
gray  scale  plot  (Figure  36)  with  Figure  39  provides  visual  (though  subjective)  evidence  of  the 
ability  of  the  ARMA  method  and  FFT  method  to  generate  equivalent  structured  scenes. 
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Figure  37.  Log-log  plot  of  two-dimensional  power  spectral  density  (PSD)  versus  spatial  frequency 
plotted  along  a  diagonal.  Curve  marked  by  an  asterisk  (•)  is  the  desired  or  "theoretical" 
PSD.  Curve  marked  b  an  X  Is  a  2-D  PSD  of  a  2-D  sheet  of  data  derived  from 
conventional  (that  is.  FFT)  simulation.  The  1024  x  1024  simulated  data  set  was  derived 
by  performing  FFT  operations  on  the  same  set  of  random  numbers  as  obtained  for 
Figure  29.  The  "theoretical"  2-D  slope  is  -2.67  and  the  correlation  length  «  84  km 
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Figure  38.  Log-log  {riot  of  one-dimensional  power  spectral  density  (PSD)  versus  spatial  frequency 
plotted  along  a  major  axis.  Curve  marked  by  an  asterisk  (•)  is  the  desired  or 
"theoretical"  1-D  PSD.  Curve  marked  by  Xs  represents  the  ID  pertodogram  derived 
from  the  rows  and  columns  of  the  2-D  FFT  simulated  data  corresponding  to  Figure  37. 
The  "theoretical"  1-D  slope  is  -1.67  and  the  correlation  length  Lq  »  84  km 
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Figure  39.  Two-dimensional  gray  scale  pictorial  representation  of  1024  x  1024  FFT  simulated  data 
(corresponds  to  Figure  37).  Hie  data  set  was  constructed  for  1^  *  84  km.  2-D  slope 
« -8/3.  spacing  ■  100  m 
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8.  COMPUTATIONAL  BURDEN 


As  mentioned  in  the  Introduction,  enormous  computer  time  is  consumed  in  constructing 
realistic  but  practical  three  dimensional  maps  of  atmospheric  temperature  or  density  fluctuations 
from  conventional  analysis  techniques.  It  is  necessary  therefore  to  explore  alternative  means  of 
constructing  structure  maps  while  reducing  the  computational  burden.  Although  the  present 
work  succeeded  in  providing  a  reliable  alternative  method,  it  did  not  realize  the  goal  of  reducing 
the  computational  burden.  This  work  does,  however,  point  to  the  possibility  of  using  a  hybrid 
FFT/AR  three-dimensional  technique  that  promises  significant  computer  savings.  This  idea  will 
be  explored  in  a  subsequent  report 

The  computer  timings  for  calculations  in  this  report  were  measured  on  the  Phillips  Laboratory 
model  210  Convex  computer.  In  solving  for  the  ARMA  filter  coefficients,  it  was  found  that  most  of 
the  time  was  spent  calculating  the  MA  terms.  For  example,  solving  for  a  3  x  3  AR  filter  took  1.2  s 
and  solving  for  a  13  x  13  AR  filter  took  3.0  s  but  solving  for  a  25  x  25  MA  filter  took  12.3  s. 
Table  1  however  shows  that  these  times  are  small  compared  to  the  time  needed  to  construct  a 
simulation  given  the  filter  coefficients.  For  example,  it  took  104  seconds  to  simulate  1300  x  1300 
data  points  using  an  order  13  x  13  AR  and  order  25  x  25  MA  filter.  In  contrast  to  the  time 
required  to  construct  the  ARMA  simulations.  Table  1  shows  that  less  than  56  seconds  were 
required  to  perform  a  conventional  2-D  FFT  simulation  corresponding  to  Figure  37  (due  to  our 
self-imposed  constraint  of  using  sections  of  the  aforementioned  ARMA  program,  approximately 
28  percent  of  the  time  is  excess  overhead).  Even  the  low  order  3  x  3  AR  and  7x7  MA  filter  took 
nearly  ten  times  as  long  to  process  as  the  conventional  method.  Thus  the  FFT  technique  proves 
much  faster  than  the  2-D  ARMA  modem  spectral  analysis  method  in  constructing  comparable  2- 
D  structure. 


Table  1.  Timing  of  Simulations 


1300  x  1300  ARMA 
(seconds) 

1024  x  1024  FFT 
(seconds) 

13  x  13  AR  (10  cycles) 

&  25  x  25  MA 

approx.  1  x  104  s 

<  56  s 

3x3  AR  (10  cycles) 

&  31x31  MA 

approx.  1  x  103  s 

<  56  s 

21  x  21  AR  (10  cycles) 
&  31  x31  MA 

approx.  3  x  104  s 

<  56  s 

3x3  AR  (20  cycles) 

&  7x  7  MA 

with  DFT  replacement 

approx.  580  s 

<  56  s 
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0.  CONCLUSIONS 


Because  geophysical  data  often  are  characterized  by  smooth  continuous  power  spectral 
density  functions,  this  report  has  explored  the  possibility  of  generating  two-dimensional  synthetic 
structure  scenes  by  passing  stochastic  data  through  an  autoregressive/moving  average  filter 
having  the  characteristics  of  the  desired  two-dimensional  PSD. 

Several  alternative  procedures  were  developed  to  construct  ARMA  filters  that  covered  a  broad 
frequency  range,  and  featured  near  circular  symmetry.  A  combination  of  quarter  plane  and  full 
plane  autoregressive  filters  cascaded  with  a  moving  average  filter  coupled  with  low  frequency 
discrete  Fourier  transform  power  replacement  proved  satisfactory  and  effective  In  generating  high 
fidelity  power  spectral  density  functions.  Simulated  In  this  way.  the  two-dimensional  power 
spectral  density  functions  of  several  scene  realizations  closely  matched  the  desired  or 
"theoretical"  PSD.  As  observed,  PSD's  having  a  small  slope,  large  correlation  length,  and  small 
data  spacing  were  more  difficult  to  model  than  PSD's  with  larger  slope,  smaller  correlation  length, 
and  wider  data  spacing.  The  analysis  Indicated  that  two-dimensional  Isotropic  horizontal  atmos¬ 
pheric  temperature  structure,  having  a  data  spacing  of  100  m  and  a  characteristically  smooth 
power  spectral  density  can  be  modeled  within  a  frequency  band  ranging  from  0.08  km1  to  5  km'1 
tty  an  autoregressive  3x3  order  process  and  a  moving  average  7x7  order  process.  DFT 
replacement  allowed  extension  of  the  lower  frequency  bound  to  0.01  km'1.  At  the  expense  of 
Increasing  computer  time,  one  may  Increase  the  frequency  range  without  DFT  replacement  and 
achieve  agreement  to  0.01  km1  by  choosing  a  13  x  13  AR  model  and  25  x  25  MA  model.  Synthetic 
scenes  generated  by  the  ARMA  process  maintained  high  fidelity  to  the  "theoretical"  PSD’s  for  two- 
dimensional  spectral  slopes  ranging  from  -8/3  to  -11/3  and  correlation  lengths  ranging  from 
32  to  84  km.  Comparison  with  a  conventional  technique  (Fast  Fourier  Transform),  showed  that 
the  ARMA  and  FFT  simulation  methods  appear  to  yield  similar  power  spectral  density  functions 
and  visual  two-dimensional  scenes. 

A  reliable  alternative  procedure  was  developed  to  produce  structure  arrays  having  circularly 
symmetric  power  spectral  densities  that  matched  desired  geophysical  specifications.  The 
procedure  accounts  for  the  isotropic  horizontal  correlation  scales.  Including  the  line-of-slght,  that 
existing  models  approximate.  The  computational  burden  of  the  ARMA  iterative  process  proved 
more  severe  than  the  conventional  FFT  technique  for  large  correlation  scales.  The  process  that 
was  eventually  chosen  led  to  run  times  10  times  the  execution  time  of  the  2-D  FFT  process. 
Despite  the  success  of  the  model  In  achieving  fidelity  to  filter  specifications,  it  is  hoped  this  report 
will  serve  as  a  detailed  study  In  2-D  ARMA  analysis  for  a  particular  geophysical  condition,  and  as 
a  caution  to  the  ability  of  method  to  achieve  computer  savings.  This  and  the  previous  work  does, 
however,  point  to  the  favorabillty  of  using  a  hybrid  FFT/AR  three-dimensional  technique  that 
promises  to  have  economical  computer  savings.  This  Idea  will  be  explored  in  a  subsequent  report 
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Appendix 

Computational  Details 


A  f  1  M  N 

x 1*  ,(m,n){cos(2^r1unfx  +  nfy )) — i  stn(2*Tj  (mfx  +  nfy ))) 


m-Mn-N 


M  N 


X  X  a,(m.  n  )(co»(2*Tj  (mfx  +  nf, )) + i  sin(2*r,  (mfx  +  nfy ))) 

m— Mn— H 

^  M  N  /  m  \2 

X  Xa,(m.n){cos[2*TI(ntfx  +  itfy)]}  +  X  ^a,(m.n){ain[2jiT,(mfx  +  nfy)]}  I 

,m»-Mn-N  )  \m— M  n— N  / 


Looking  at  the  numerators,  Nj  and  N4  of  — —  and  — — ,  we  observe, 

*ai  aA4 

M  N 

,(m,n){cos[2nT1(mfx  +  nfy }]  + 1  sfn[2xT,(mfx  +  nfy)]} 

m«0n>0 

0  0,  , 

N4  *  £  ^a4(m,n){cosf 2rtTl(mfx  +  nfy )]  + 1  sln[2itT1  (mfx  +  nfy)]| 

m=-M n«-N 

but  a,(m.n)  =  a4(-m,-n) 
so. 

N,  =  tia  4(m,n){cos[2RTi(-mfx  -nfy)]+lsln[27<T1(-infx  -nfy)]} 
n>0n>0 

so  that. 

M  N 

N,  +  N4=2XXa  ,(m.n)cos[2jtT,(mfx  +  nfy )] 

m«0n«0 
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Looking  at  the  denominators  D,  and  D4  of  — —  and  — ,  we  observe. 

AA1  AA4 


D,  = 


D4  = 


M  N 


m«0n«0 
0  0 


M  N 

^^a,(m.n)cos[2?tT,(infx  +nfy)]^  +  j  ^£a,(m.n)sin[2tfr,(mfx  +nfy)] 

n*0n*0 


lm>-Mn>- N 

but.  a,(m.n)  =  a4(-m.-n).  so  that. 


00 

£  £a4(m.n)cos[2nT,(mfx+nfy)]}  +j  £  ^a4{m.n)sln[2flT,(mfx +  nfy)] 


(  M  N 


Da  = 


Ym=0n=0 


(  M  N 


£  ^a1(m.n)cos(2«T1(-mfx  -  nfy ))  +  IS  a1(m.n)stn(2jcT,(-mfx  -nfy)) 


\m*0n*0 


so  that.  D,  =  D4 
consequently. 

M  N 

.  .  ’ll*  1(m,n)cos[2jtT1(mfx  +  nfy)j 

— .  +  — L.  - _ maOn»Q _ _ _ 

M  N 

IS*  j(m,n)sln[27tT,(mfx  +  nfy)] 

m»0n«0 


'Al 


*A4 


M  N  V 

SI-  ,(m.n)cos[2nT,(mfx  +  nfy)H 
.m=0n*0  J 


•H 


Likewise  looking  at  the  numerators.  N2  and  N3  of  — - —  and  — - — .  we  observe. 

Aa2  Aa3 

ON,  . 

N2=  £  5rfaa(m,n){cos[2sT|(mfx  +  nfy )]  + 1  sta[2ifr1(mfx  +  nfy)]j 

m=-M  n=0 

but  a2(-m,n)  =  a3(m.-n) 

So. 


M  N  . 

n2  =  XXa2(  -m,n)jcos[2sTj(-mfx  +  nfy)]+lsin[2rtTi(-mfx  +nfy) 

m=0n=0 

M  N  f  , 

N3  =  IS*  3(m.-n){cos[2jtT,(mfx  -nfy)]  +  isln[2jiT,(mfx  -nfy)]} 

m=0n=0 

but  a3(m,-n)  =  a2(-m,n)  =  a,(m,n) 
so  that. 

M  N 

N2  +  N3  =  2  IS*  i(m,n)cos[2jtTj{mfx  -  nfy)] 

m=0n=0 
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Looking  at  the  denominators  D2  and  D3  of  — and  — ,  we  observe. 

“A2  “A3 


Da  = 

Dj  = 


ON  la  r  0  N  Y 

£a2(m.n)cos[2rfr,(mfx+nfy)H  +j  X  ]£a2(m.n)sln[2rfr,(mfx  +  nfy)H 
i--M  n-0  J  (.m— M  n-0  J 

'  M  N  l2  I"  M  N  l2 

X  Xa2*“m,n^^2rfri*~infx  +nfy*][  +  |  X  Xa2(-m*n)aln[2,n'i('mf*  +  nfy)]f 

m-0  n-0  J  1,  m-0  n-0  J 

MO  |3  f  M  0  l2 

D3=  £  ^a3(m, n)cos[2rfT,  (mfx  +  nfy  )H  +|X  Xa3fm,n*sin[2,rr^inf* +nfyj][ 
m-0  n— N  J  (.m-On— N  J 

[  M  N  l2  f  M  N  l2 

D3  =  ]X  Xa3(m.-n)c°s[2,fr1(infx-nfy)]  +  j  £  ^a3(m.-n)sin[2flT,(mfx -nfy)H 
Im-O  n-0  J  [m-0  n-0  J 


but  a2(-m,n)  =  a3(m,-n)  =  a^m.n) 

so  that.  D2  =  D3 

consequently, 


M  N 


1  t  1 
AA2  AA3 


Finally 


2^T  ]£a1(m.n)cos[2i«T,(mfx-nfy)] 

_ m-0  n-0 _ _ _ 


M  N  p  M  N  p 

(m.n)co^2rfT,(mfx  -nfy)]>  +1  ^  ^a  lm.n)sin[2nTIlmfx  -nfy)]> 

m-On-O  lm<0n<0 


M  N 


2Xy.alcosf2*Tlfrnf»4nfy^ 

n>OnC _ _ _ 


M  N 


2yy.a,cosf2iff,(irfx-nfj| 
artofl _ 


fMN  '  *  M  N  '  *  ’  M  N  3  f|<  N  '  • 

| XXH2 *TiM,+n9]  +  XXjsln[2j<rj(in4-m^)j  -  Y.YcQ8f2i<rltnrf,-rfv)]  +|^^sln[2*T1(ni4-nfy)j 

UnOnO  UDOnO  IptOmO  UdOdO 
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Ignoring  the  sine  terms  for  fx  -*  0,  and  fy  -» 0 : 


£  XX-  i(m,n)cos[2xT,(mfx  +nfy)j 

£-*o  na0n>0 


M  IS 

XX-  ,(m,n)[cos(2xmfxT1)cos(2xnfyT1)  -  flta(2xmfxT,)sin(2imfyT,)] 


m»0  n»0 


m  w 

XX*  ,  (m.n)[co8(2xmfxT,  )cos(2xnfyT, )] 


m»0  n-0 


where,  again  we  have  ignored  sine  terms.  This  Implies  that  as  fx  ->  0. and  fy  -+  0,  and  checking 

1  4 

the  definition  of  al(m.n),  we  see  as  before,  the  following  approximation:  — — — — -  -  ■  It  Is 

unclear  why  the  sum  of  the  quarter  plane  filters  performs  better  than  the  average. 

ANALYSIS  OF  LOW  FREQUENCY  DFT  REPLACEMENT 

The  following  describes  the  analysis  and  procedure  for  subtracting  the  low 
frequency  components  of  the  2-D  ARMA  simulated  scene. 

Assume  that  G(x,y)  is  a  two-dimensional  sheet  of  simulated  data.  We  wish  to  find  the  discrete 
Fourier  coefficients  Alaij.o^)  at  frequencies  o>i  and  »2  where. 


but  not  to,  =  co2  =  0  or  |<o,|  =  jcc^l  =  x.  Then. 

Afw,.a>a; = ^■X2G(x,y)e’tahXe'tosy 


A(-o),.-a)2)  =  ^-J^Glx.yle^V*** 


=s  A*(C0,.C02) 
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So,  where  9)  signifies  the  real  components,  and  3  signifies  the  Imaginary  components: 

StAjoOi .  eoj)  -  SKaJ-oj,  ,  -coj) 

*  ■  ^^^CCx.y^eos^tOjxJa^to^y)- sinfeojxJsinloajy)] 

x  y 

,  and. 

SAtOpOOs)  =  -SAl-c^.-tOa) 

=  T?r^5)GCx.y)[cO6(<0jx)sin(tD2y)+  sinltojxk^sloojy)] 
x  y 

Now  the  contribution  of  A(a>j  .o^)  at  point  x.y  is: 

C(x.y)  =  9?[A(©,.o^)][cos(co,x)cos(ci>2y)  -  sinftOjXjainftOay)] 

+3[A(a>1,tfla)][co9((o1x)sln((i>2y)  +  aln((OiX)co8((o2y)] 

and  the  contribution  of  AC-ai!  .-aij)  at  point  x,y  is: 

=  9?[A(-coj  .-to2)][cos(coix}cos(Q)2y)  -  8ln((o,x)sln(co2y)] 

+  3[  A(-ca, .  -0)3  )][-cos((o,x  )sin(w2y )  -  slnfb^xjcostci^ )] 

and  since  Afto^cOj)  =  A*(-tolf-co2).  the  contribution  of  A(<d,,(D2)  +  A(-0)lt-0)2)  is : 

CC(x.y)  =  29?[A(io,,co2)Jcos(io1x)cos((o2y)-  stnfo^xJslnCa^y)] 

+23jA(a>, ,  <o2  )|cos(co,x)  stn(©2y)  +  stnlwjxlcosCa^y )] 

Taking  the  subset  of  frequencies: 

-110^-jstOj  and  0  £  ©2  £  1  but  not  ©j  =  ©3  =  0 

(we  call  this  a  12  x  12  or  L  x  L  subset  of  frequencies)  we  subtract  CC(x,y)  at  these  frequencies 
from  G(x,y)  and  add  instead  the  DD(x,y)  values  described  below. 
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The  following  describes  the  analysis  and  procedure  for  replacing  the  low 
frequency  components  to  the  simulated  scene. 


We  set  up  an  array  of  2L  x  2L  pseudo- random  numbers  efs.t)  drawn  from  a  Gaussian 
probability  distribution  with  standard  deviation  *  1  and  then  calculate  the  2-D  FFT  of  e(s,t),  Gfa^, 
0)4).  where. 


a>3  = 

a»4  = 


2xT 

2L 

2xT 

2L 


-LSTSL. 


Invoking  Pascal's  theorem,  the  expected  value  of  the  square  of  the  amplitudes  of  Gta^.co^  is  (2LJ2. 
But  G(cf>3.a>4)  has  the  proper  functional  form  of  the  probability  density  function  for  the  Fourier 
transform  of  the  data  we  wish  to  simulate.  The  Idea  Is  to  take  a  value  of  and  multiply  It 

by  a  value  Mfcoj.oj)  to  get  the  expected  value  of  Afcoj.a^)  which  will  have  the  power  spectral 
density  we  wish  to  simulate.  This  means: 

|GG(o>3,(o4)|2[M(o>,.(a2)]2  =  w^ere  Ffwj.ooj)  is  the  power  spectral  density  we  are  trying  to 

simulate  and  Ax  •  100  m  Is  the  spacing  of  the  data  we  are  simulating.  Thus. 


or. 


M(<Ui.<o2) 


jF((ot.(oa) 
2x AxxLxN 


and. 


Dfwj.coa)  =  GGf^-.^-l'Mfo)!.^) 

zL  zL 

Then,  as  above,  the  contribution  at  point  (x,y)  of  DfCDptDa) +  D(-C0|,-(D2)  is  : 
DD(x.y)  =  29t[D(a>1,to2)Jcos(co1x)cos(to2y)  -  slnf^xlslnfa^y)] 

+23[D(©1  .cua  )][cos(a),x)sin((o2y ) + slnfa^xlcosCa^y)] 


» 
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